library(TSA)
library(tseries)
library(astsa)
library(imputeTS)
library(tsoutliers)
library(xts)
Here we interpolate our missing data with a linear model.
terror2 <- read.csv("input/og_num_casualities_greater_than_10.csv")
terror3 <- na.interpolation(terror2$num.attacks.with.kill.thresh, option="linear")
plot(as.xts(ts(terror3, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (w/ Linear Imputed Data)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", ylim=c(0, 225), col="red")
pdf("image/og_ts.pdf")
lines(as.xts(ts(terror2$num.attacks.with.kill.thresh, frequency = 12, start=1970)), col="black", lwd=1.2)
dev.off()
quartz_off_screen
2
outlier_terror3 <- tso(ts(terror3), types = c("TC", "AO", "IO"))
stopped when 'maxit.oloop = 4' was reached
plot(outlier_terror3)
#plot outlier effects
pdf("image/outlier_effects.pdf")
plot(as.xts(ts(outlier_terror3$effects, frequency = 12, start=1970)), main = "Outlier Effects", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="red")
dev.off()
quartz_off_screen
2
#Plot outlier time series
xts.terror3 <- as.xts(ts(terror3, frequency = 12, start=1970))
plot(as.xts(ts(terror3, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Outliers Removed)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="lightgray", ylim=c(0, 225))
lines(as.xts(ts(outlier_terror3$yadj, frequency = 12, start=1970)), col="blue")
points(xts.terror3[427], col="red",pch=19, cex=1)
points(xts.terror3[516], col="red",pch=19, cex=1)
points(xts.terror3[521], col="red",pch=19, cex=1)
points(xts.terror3[523], col="red",pch=19, cex=1)
points(xts.terror3[547], col="red",pch=19, cex=1)
pdf("image/outlier_comparison.pdf")
points(xts.terror3[556], col="red",pch=19, cex=1)
dev.off()
quartz_off_screen
2
terror4 <- outlier_terror3$yadj
#terror3 <- na.kalman(terror2$num.attacks, model="auto.arima")
cuttoff.index <- length(terror4) - 24 #floor(0.1 * length(terror3))
cuttoff.index2 <- length(terror4) - 12
terror4.valid <- terror4[(cuttoff.index+1) :cuttoff.index2]
terror4.testing <- terror4[(cuttoff.index2 + 1): length(terror4)]
terror4 <- terror4[1: cuttoff.index]
#plot(as.xts(ts(terror4, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Training Set)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
#plot(as.xts(ts(terror4.valid, frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Validation Set)", major.format = "%Y-%m", grid.col="white", lwd=1)
#log_terror4 <- log(outlier_terror3$yadj)
adf.test(terror4, k=1)
Augmented Dickey-Fuller Test
data: terror4
Dickey-Fuller = -2.4189, Lag order = 1, p-value = 0.401
alternative hypothesis: stationary
adf.test(diff(terror4), k=1)
p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: diff(terror4)
Dickey-Fuller = -24.444, Lag order = 1, p-value = 0.01
alternative hypothesis: stationary
adf.test(diff(diff(terror4)), k=1)
p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: diff(diff(terror4))
Dickey-Fuller = -31.986, Lag order = 1, p-value = 0.01
alternative hypothesis: stationary
pdf("image/first_diff.pdf")
plot(as.xts(ts(diff(terror4), frequency = 12, start=1970)), main = "Number of Terrorist Attacks (First Diff)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
dev.off()
null device
1
pdf("image/second_diff.pdf")
plot(as.xts(ts(diff(diff(terror4)), frequency = 12, start=1970)), main = "Number of Terrorist Attacks (Second Diff)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years")
dev.off()
null device
1
#ts.plot(diff(terror4))
#ts.plot(diff(diff(terror4)))
pdf("image/acf_og.pdf")
acf(terror4, main="ACF of Training Data")
dev.off()
null device
1
pdf("image/pacf_og.pdf")
pacf(terror4, main="PACF of Training Data")
dev.off()
null device
1
pdf("image/acf_first_diff.pdf")
acf(diff(terror4), main="ACF of First Diff Training Data")
dev.off()
null device
1
pdf("image/pacf_first_diff.pdf")
pacf(diff(terror4), main="PACF of First Diff Training Data")
dev.off()
null device
1
pdf("image/acf_second_diff.pdf")
acf(diff(terror4), main="ACF of Second Diff Training Data")
dev.off()
null device
1
pdf("image/pacf_second_diff.pdf")
pacf(diff(terror4), main="PACF of Second Diff Training Data")
dev.off()
null device
1
m = floor(sqrt(length(diff(terror4))))
#pdf("image/raw_periodogram.pdf")
spec.pgram(diff(terror4), log="no", main="Raw Periodogram", cex.main=1.5)
#dev.off()
#pdf("image/smooth_tapered_periodogram.pdf")
spec.pgram(diff(terror4), kernel('daniell', m), log="no", taper=0.1, main="Smoothed and Tapered Periodogram")
#dev.off()
#mvspec(diff(log_terror4), kernel('daniell', m), log="no")
eacf(diff(terror4))
AR/MA
0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o x x o o o o o o o o o x
1 x o o x o o o o o o o o o o
2 x x x x o o o o o o o o x o
3 o x x x o o o o o o o o x o
4 x o o o x o x o o o o o o o
5 x x o o o o o o o o o o o o
6 x x x o o o o o o o o o o o
7 x x x x o o o o o o o o o o
eacf(diff(diff(terror4)))
AR/MA
0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x x x x o o o o o o o o o x
1 x x x x o o o o o o o o o x
2 x x x x o x o o o o o o o o
3 x x o o o o o o o o o o o o
4 x x o o o o x x o o o o o o
5 x x o o o o o o o o o o o o
6 x o o x o o o o o o o o o o
7 x x o x x o o o o o o o o o
#sarima(terror4, 0, 1, 1)
sarima(terror4, 0, 1, 1, 1, 0, 1, 4)
initial value 2.346947
iter 2 value 2.220389
iter 3 value 2.211951
iter 4 value 2.193176
iter 5 value 2.187658
iter 6 value 2.185821
iter 7 value 2.185525
iter 8 value 2.185512
iter 9 value 2.185508
iter 10 value 2.185484
iter 11 value 2.185474
iter 12 value 2.185467
iter 13 value 2.185459
iter 14 value 2.185419
iter 15 value 2.185268
iter 16 value 2.184298
iter 17 value 2.183991
iter 18 value 2.182728
iter 19 value 2.181869
iter 20 value 2.181763
iter 21 value 2.181685
iter 22 value 2.181580
iter 23 value 2.181485
iter 24 value 2.181480
iter 25 value 2.181480
iter 25 value 2.181480
iter 25 value 2.181480
final value 2.181480
converged
initial value 2.178465
iter 2 value 2.178461
iter 3 value 2.178459
iter 4 value 2.178459
iter 5 value 2.178458
iter 6 value 2.178457
iter 7 value 2.178457
iter 8 value 2.178457
iter 8 value 2.178457
iter 8 value 2.178457
final value 2.178457
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))
Coefficients:
ma1 sar1 sma1 constant
-0.6128 -0.6221 0.7609 0.2964
s.e. 0.0351 0.1214 0.1003 0.1602
sigma^2 estimated as 77.91: log likelihood = -1939, aic = 3887.99
$degrees_of_freedom
[1] 535
$ttable
Estimate SE t.value p.value
ma1 -0.6128 0.0351 -17.4617 0.0000
sar1 -0.6221 0.1214 -5.1236 0.0000
sma1 0.7609 0.1003 7.5870 0.0000
constant 0.2964 0.1602 1.8496 0.0649
$AIC
[1] 5.370387
$AICc
[1] 5.374299
$BIC
[1] 4.402176
sarima(terror4, 0, 1, 1, 1, 1, 1, 4)
initial value 2.599103
iter 2 value 2.332712
iter 3 value 2.270429
iter 4 value 2.237537
iter 5 value 2.222895
iter 6 value 2.215065
iter 7 value 2.205793
iter 8 value 2.204129
iter 9 value 2.203212
iter 10 value 2.203071
iter 11 value 2.202956
iter 12 value 2.202950
iter 13 value 2.202949
iter 13 value 2.202949
iter 13 value 2.202949
final value 2.202949
converged
initial value 2.204741
iter 2 value 2.203758
iter 3 value 2.203347
iter 4 value 2.203343
iter 5 value 2.203343
iter 6 value 2.203342
iter 6 value 2.203342
iter 6 value 2.203342
final value 2.203342
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ma1 sar1 sma1
-0.6121 0.1591 -1.0000
s.e. 0.0375 0.0434 0.0362
sigma^2 estimated as 79.13: log likelihood = -1937.92, aic = 3883.84
$degrees_of_freedom
[1] 532
$ttable
Estimate SE t.value p.value
ma1 -0.6121 0.0375 -16.3203 0e+00
sar1 0.1591 0.0434 3.6675 3e-04
sma1 -1.0000 0.0362 -27.6028 0e+00
$AIC
[1] 5.382182
$AICc
[1] 5.386024
$BIC
[1] 4.406024
sarima(terror4, 0, 1, 1, 1, 1, 2, 4)
initial value 2.599103
iter 2 value 2.300969
iter 3 value 2.258534
iter 4 value 2.236661
iter 5 value 2.205300
iter 6 value 2.201253
iter 7 value 2.198386
iter 8 value 2.197243
iter 9 value 2.196562
iter 10 value 2.196267
iter 11 value 2.196213
iter 12 value 2.196133
iter 13 value 2.195850
iter 14 value 2.194818
iter 15 value 2.194074
iter 16 value 2.193350
iter 17 value 2.192729
iter 18 value 2.192231
iter 19 value 2.192117
iter 20 value 2.192100
iter 21 value 2.192094
iter 22 value 2.192032
iter 23 value 2.191993
iter 24 value 2.191956
iter 25 value 2.191945
iter 26 value 2.191944
iter 26 value 2.191944
iter 26 value 2.191944
final value 2.191944
converged
initial value 2.199972
iter 2 value 2.199945
iter 3 value 2.199936
iter 4 value 2.199929
iter 5 value 2.199892
iter 6 value 2.199876
iter 7 value 2.199872
iter 8 value 2.199872
iter 8 value 2.199872
iter 8 value 2.199872
final value 2.199872
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ma1 sar1 sma1 sma2
-0.6168 -0.6197 -0.2244 -0.7477
s.e. 0.0381 0.1241 0.1090 0.0999
sigma^2 estimated as 79.3: log likelihood = -1936.06, aic = 3882.13
$degrees_of_freedom
[1] 531
$ttable
Estimate SE t.value p.value
ma1 -0.6168 0.0381 -16.2046 0.00
sar1 -0.6197 0.1241 -4.9951 0.00
sma1 -0.2244 0.1090 -2.0593 0.04
sma2 -0.7477 0.0999 -7.4861 0.00
$AIC
[1] 5.388087
$AICc
[1] 5.391999
$BIC
[1] 4.419876
sarima(terror4, 1, 1, 2)
initial value 2.344191
iter 2 value 2.241534
iter 3 value 2.195978
iter 4 value 2.194817
iter 5 value 2.194230
iter 6 value 2.194208
iter 7 value 2.194208
iter 8 value 2.194207
iter 9 value 2.194207
iter 10 value 2.194205
iter 11 value 2.194205
iter 12 value 2.194204
iter 13 value 2.194204
iter 14 value 2.194204
iter 15 value 2.194201
iter 16 value 2.194195
iter 17 value 2.194183
iter 18 value 2.194180
iter 19 value 2.194178
iter 20 value 2.194178
iter 21 value 2.194178
iter 21 value 2.194178
iter 21 value 2.194178
final value 2.194178
converged
initial value 2.193663
iter 2 value 2.193663
iter 3 value 2.193663
iter 3 value 2.193663
iter 3 value 2.193663
final value 2.193663
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))
Coefficients:
ar1 ma1 ma2 constant
-0.0950 -0.4967 -0.0688 0.2960
s.e. 0.7203 0.7171 0.4282 0.1536
sigma^2 estimated as 80.36: log likelihood = -1947.19, aic = 3904.38
$degrees_of_freedom
[1] 535
$ttable
Estimate SE t.value p.value
ar1 -0.0950 0.7203 -0.1320 0.8951
ma1 -0.4967 0.7171 -0.6926 0.4888
ma2 -0.0688 0.4282 -0.1608 0.8723
constant 0.2960 0.1536 1.9265 0.0546
$AIC
[1] 5.401315
$AICc
[1] 5.405227
$BIC
[1] 4.433105
sarima(terror4, 1, 1, 2, 1, 0, 1, 4)
initial value 2.347868
iter 2 value 2.258620
iter 3 value 2.188794
iter 4 value 2.186745
iter 5 value 2.185962
iter 6 value 2.185836
iter 7 value 2.185771
iter 8 value 2.185767
iter 9 value 2.185742
iter 10 value 2.185679
iter 11 value 2.185479
iter 12 value 2.185171
iter 13 value 2.184101
iter 14 value 2.183062
iter 15 value 2.182366
iter 16 value 2.182283
iter 17 value 2.182258
iter 18 value 2.182224
iter 19 value 2.182152
iter 20 value 2.182073
iter 21 value 2.182032
iter 22 value 2.181984
iter 23 value 2.181979
iter 24 value 2.181979
iter 24 value 2.181979
iter 24 value 2.181979
final value 2.181979
converged
initial value 2.178050
iter 2 value 2.178045
iter 3 value 2.178044
iter 4 value 2.178044
iter 5 value 2.178044
iter 6 value 2.178043
iter 7 value 2.178043
iter 8 value 2.178043
iter 8 value 2.178043
iter 8 value 2.178043
final value 2.178043
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))
Coefficients:
ar1 ma1 ma2 sar1 sma1 constant
0.1981 -0.7935 0.0937 -0.6070 0.7515 0.2966
s.e. 0.8900 0.8938 0.5535 0.1345 0.1104 0.1555
sigma^2 estimated as 77.85: log likelihood = -1938.77, aic = 3891.55
$degrees_of_freedom
[1] 533
$ttable
Estimate SE t.value p.value
ar1 0.1981 0.8900 0.2226 0.8239
ma1 -0.7935 0.8938 -0.8878 0.3750
ma2 0.0937 0.5535 0.1692 0.8657
sar1 -0.6070 0.1345 -4.5117 0.0000
sma1 0.7515 0.1104 6.8085 0.0000
constant 0.2966 0.1555 1.9074 0.0570
$AIC
[1] 5.376974
$AICc
[1] 5.381068
$BIC
[1] 4.424658
sarima(terror4, 1, 1, 2, 1, 1, 1, 4)
initial value 2.600040
iter 2 value 2.356718
iter 3 value 2.277980
iter 4 value 2.233407
iter 5 value 2.223552
iter 6 value 2.212201
iter 7 value 2.205827
iter 8 value 2.203297
iter 9 value 2.201565
iter 10 value 2.200773
iter 11 value 2.200719
iter 12 value 2.200705
iter 13 value 2.200661
iter 14 value 2.200544
iter 15 value 2.200182
iter 16 value 2.200058
iter 17 value 2.199536
iter 18 value 2.198925
iter 19 value 2.198604
iter 20 value 2.198391
iter 21 value 2.198352
iter 22 value 2.198279
iter 23 value 2.198233
iter 24 value 2.197691
iter 25 value 2.195332
iter 26 value 2.195248
iter 27 value 2.193232
iter 28 value 2.191968
iter 29 value 2.190157
iter 30 value 2.190012
iter 31 value 2.189751
iter 32 value 2.189563
iter 33 value 2.189216
iter 34 value 2.188867
iter 35 value 2.188828
iter 36 value 2.188810
iter 37 value 2.188802
iter 37 value 2.188802
iter 37 value 2.188802
final value 2.188802
converged
initial value 2.202122
iter 2 value 2.202104
iter 3 value 2.202047
iter 4 value 2.202036
iter 5 value 2.202000
iter 6 value 2.201960
iter 7 value 2.201834
iter 8 value 2.201715
iter 9 value 2.201577
iter 10 value 2.201560
iter 11 value 2.201548
iter 12 value 2.201543
iter 12 value 2.201543
final value 2.201543
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ar1 ma1 ma2 sar1 sma1
0.4461 -1.0292 0.2206 0.1873 -0.9933
s.e. 0.3240 0.3308 0.2148 0.0489 0.0407
sigma^2 estimated as 79.29: log likelihood = -1936.96, aic = 3885.92
$degrees_of_freedom
[1] 530
$ttable
Estimate SE t.value p.value
ar1 0.4461 0.3240 1.3767 0.1692
ma1 -1.0292 0.3308 -3.1116 0.0020
ma2 0.2206 0.2148 1.0270 0.3049
sar1 0.1873 0.0489 3.8268 0.0001
sma1 -0.9933 0.0407 -24.3931 0.0000
$AIC
[1] 5.39168
$AICc
[1] 5.395676
$BIC
[1] 4.431417
sarima(terror4, 1, 1, 2, 1, 1, 2, 4)
initial value 2.600040
iter 2 value 2.325435
iter 3 value 2.237416
iter 4 value 2.214251
iter 5 value 2.202589
iter 6 value 2.198408
iter 7 value 2.198339
iter 8 value 2.197318
iter 9 value 2.196580
iter 10 value 2.196494
iter 11 value 2.196379
iter 12 value 2.196322
iter 13 value 2.195990
iter 14 value 2.195374
iter 15 value 2.195370
iter 16 value 2.195119
iter 17 value 2.194986
iter 18 value 2.194830
iter 19 value 2.194773
iter 20 value 2.194505
iter 21 value 2.194160
iter 22 value 2.193795
iter 23 value 2.193655
iter 24 value 2.193455
iter 25 value 2.193431
iter 26 value 2.193425
iter 27 value 2.193424
iter 27 value 2.193424
final value 2.193424
converged
initial value 2.199628
iter 2 value 2.199619
iter 3 value 2.199609
iter 4 value 2.199604
iter 5 value 2.199600
iter 6 value 2.199580
iter 7 value 2.199565
iter 8 value 2.199563
iter 9 value 2.199562
iter 10 value 2.199561
iter 11 value 2.199557
iter 12 value 2.199550
iter 13 value 2.199534
iter 14 value 2.199505
iter 15 value 2.199481
iter 16 value 2.199471
iter 17 value 2.199467
iter 18 value 2.199467
iter 19 value 2.199467
iter 20 value 2.199464
iter 21 value 2.199458
iter 22 value 2.199449
iter 23 value 2.199441
iter 24 value 2.199438
iter 25 value 2.199436
iter 26 value 2.199433
iter 27 value 2.199426
iter 28 value 2.199422
iter 29 value 2.199420
iter 30 value 2.199419
iter 31 value 2.199417
iter 32 value 2.199414
iter 33 value 2.199410
iter 34 value 2.199407
iter 35 value 2.199407
iter 36 value 2.199407
iter 37 value 2.199406
iter 38 value 2.199406
iter 39 value 2.199404
iter 40 value 2.199400
iter 40 value 2.199400
final value 2.199400
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ar1 ma1 ma2 sar1 sma1 sma2
0.0787 -0.6766 0.0177 -0.6062 -0.2298 -0.7358
s.e. 0.9149 0.9129 0.5568 0.1362 0.1180 0.1096
sigma^2 estimated as 79.34: log likelihood = -1935.81, aic = 3885.62
$degrees_of_freedom
[1] 529
$ttable
Estimate SE t.value p.value
ar1 0.0787 0.9149 0.0860 0.9315
ma1 -0.6766 0.9129 -0.7412 0.4589
ma2 0.0177 0.5568 0.0318 0.9746
sar1 -0.6062 0.1362 -4.4501 0.0000
sma1 -0.2298 0.1180 -1.9474 0.0520
sma2 -0.7358 0.1096 -6.7148 0.0000
$AIC
[1] 5.395983
$AICc
[1] 5.400077
$BIC
[1] 4.443667
sarima(terror4, 1, 1, 1)
initial value 2.344191
iter 2 value 2.241424
iter 3 value 2.210360
iter 4 value 2.203987
iter 5 value 2.203217
iter 6 value 2.196675
iter 7 value 2.194860
iter 8 value 2.194317
iter 9 value 2.194264
iter 10 value 2.194223
iter 11 value 2.194209
iter 12 value 2.194204
iter 13 value 2.194204
iter 13 value 2.194204
iter 13 value 2.194204
final value 2.194204
converged
initial value 2.193680
iter 2 value 2.193680
iter 3 value 2.193679
iter 3 value 2.193679
iter 3 value 2.193679
final value 2.193679
converged
$fit
Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))
Coefficients:
ar1 ma1 constant
0.0175 -0.6101 0.2960
s.e. 0.0678 0.0517 0.1537
sigma^2 estimated as 80.36: log likelihood = -1947.2, aic = 3902.4
$degrees_of_freedom
[1] 536
$ttable
Estimate SE t.value p.value
ar1 0.0175 0.0678 0.2580 0.7965
ma1 -0.6101 0.0517 -11.7941 0.0000
constant 0.2960 0.1537 1.9258 0.0547
$AIC
[1] 5.397646
$AICc
[1] 5.401488
$BIC
[1] 4.421488
pdf("image/best_model.pdf")
sarima(terror4, 0, 1, 1, 1, 1, 1, 4)
dev.off()
#eacf(diff(diff(log_terror4)))
#terror5 <- terror4
#total_error <- 0
#for (i in 1: (length(terror4.valid) - 11))
#{
# actual <- terror4.valid[i : i + 11]
# predicted <- sarima.for(terror5, 12, 0, 1, 1, 0, 0, 0, 0) $pred
# total_error <- total_error + sum((actual - predicted)^2)
# terror5 <- c(terror5, terror4.valid[i])
#}
predicted <- sarima.for(terror4, 12, 0, 1, 1, 0, 0, 0, 0)$pred
mse <- sum((predicted - terror4.valid)^2)
mse
[1] 4248.281
val <- sarima.for(c(terror4, terror4.valid), 12, 0, 1, 1, 1, 1, 1, 4)
pred <-val$pred
err <-val$se
total <- c(terror4, terror4.valid, terror4.testing)
par(cex.main = 2)
plot(as.xts(ts(total, frequency = 12, start=1970))[492:length(total)], main = "Number of Terrorist Attacks (Prediction)", major.format = "%Y-%m", grid.col="white", lwd=1, major.ticks = "years", col="lightgray", pch="1", ylim=c(0, 225))
points(as.xts(ts(total, frequency = 12, start=1970)),col="lightgray",pch="o")
lines(as.xts(ts(c(terror4, terror4.valid), frequency = 12, start=1970)),col="black")
points(as.xts(ts(c(terror4, terror4.valid), frequency = 12, start=1970)),col="black",pch="o")
lines(as.xts(ts(pred, frequency = 12, start=2016)),col="blue")
lines(as.xts(ts(pred + err, frequency = 12, start=2016)),col="blue", lty="dashed")
lines(as.xts(ts(pred - err, frequency = 12, start=2016)),col="blue", lty="dashed")
lines(as.xts(ts(pred + 2*err, frequency = 12, start=2016)),col="blue", lty="dotted")
pdf("image/prediction_on_testing.pdf")
lines(as.xts(ts(pred - 2*err, frequency = 12, start=2016)),col="blue", lty="dotted")
dev.off()
quartz_off_screen
2
mse <- sum((pred - terror4.testing)^2)
mse
[1] 1784.668